remove(list=ls())

setwd("/Users/cb2257/Desktop/APSR Replication/")

library("lme4")
library("bayestestR")
library("tidyverse")
library("dplyr")
library("leaflet")
library("haven")
library("leaflet")
library("ggplot2")
library("mosaic")
library("estimatr")

data <- read_dta("Data/Survey/survey_final.dta")
estsample <- data[which(data$estsample==1),]

estsample <- estsample %>% mutate(std_age = scale(age))
estsample <- estsample %>% mutate(std_ideology = scale(ideology))

varlist <- c("republican", "democrat", "woman", "highschool", "college", "std_age", "std_ideology", "race_white", "race_black", "race_latinx", "employed", "native", "religiosity", "empathy", "pol_interest")

#Perform Equivalence Tests -- Figure A-2 in Main Appendix

est_equiv <- function(varname) {
require(parameters)
require(stats)
  
#Define Empirical Model
f <- paste0("hurricane_index1 ~ factor(start_date) + factor(statefp) + ", varname)

#Estimate Model With Weights
m.wt <- lm(f, data = estsample, weights = eweight)

#Conduct Equivalence Test
equiv.wt <- equivalence_test(m.wt, rule = "classic", ci=.975, range(c(-.15,.15)))

##Extract Model Information
df.wt <- data.frame(equiv.wt)
df.wt$sample <- " "
df.wt$estimate <- coef(m.wt)

#subset intercept since not of interest
df.out <- subset(df.wt, Parameter!=c("(Intercept)", "factor(start_date)2022-08-22" , "factor(start_date)2022-08-23", "factor(start_date)2022-08-24", "factor(start_date)2022-08-26", "factor(start_date)2022-08-27", "factor(start_date)2022-08-29", "factor(start_date)2022-09-06", "factor(start_date)2022-09-19", "factor(start_date)2022-09-26", "factor(start_date)2022-09-27", "factor(start_date)2022-09-28", "factor(start_date)2022-09-29", "factor(start_date)2022-10-03", "factor(start_date)2022-10-04", "factor(start_date)2022-10-07", "factor(start_date)2022-10-12", "factor(start_date)2022-10-21", "factor(start_date)2022-10-27"))
return(df.out)
}
#loop over variable list
equiv.out <- lapply(varlist, est_equiv)
equiv.out <- do.call("rbind", equiv.out)
equiv.out

republican <- subset(equiv.out,Parameter=="republican")
democrat <- subset(equiv.out,Parameter=="democrat")
highschool <- subset(equiv.out,Parameter=="highschool")
college <- subset(equiv.out,Parameter=="college")
woman <- subset(equiv.out,Parameter=="woman")
age <- subset(equiv.out,Parameter=="std_age")
ideology <- subset(equiv.out,Parameter=="std_ideology")
white <- subset(equiv.out,Parameter=="race_white")
black <- subset(equiv.out,Parameter=="race_black")
latinx <- subset(equiv.out,Parameter=="race_latinx")
employed <- subset(equiv.out,Parameter=="employed")
native <- subset(equiv.out,Parameter=="native")
religiosity <- subset(equiv.out,Parameter=="religiosity")
empathy <- subset(equiv.out,Parameter=="empathy")
interest <- subset(equiv.out,Parameter=="pol_interest")

equiv.out<-rbind(republican, democrat, highschool, college, woman, age, ideology, white, black, latinx, employed, native, religiosity, empathy, interest)

equiv.out$Parameter  = factor(equiv.out$Parameter, levels=c("pol_interest", "empathy", "religiosity", "native", "employed", "race_latinx", "race_black", "race_white", "std_ideology", "std_age", "college", "highschool", "woman", "democrat", "republican"))

equiv.out$Estimate <- equiv.out$sample
equiv.out$Hurricane <- equiv.out$estimate 

renames <- c("republican" = "Republican", "democrat" = "Democrat", "highschool" = "High School Graduate", "college" = "College Graduate", "woman" = "Woman", "std_age" = "Age", "std_ideology" = "Political Ideology", "race_white" = "White", "race_black" = "Black", "race_latinx" = "Latinx", "employed" = "Employed", "native" = "Native Born", "religiosity" = "Religiosity", "empathy" = "Empathy", "pol_interest" = "Political Interest")
equiv.out <- equiv.out %>% mutate(Parameter = str_replace_all(Parameter, renames))

#create balance plot
rope_low <- equiv.out$ROPE_low[[1]]
rope_high <- equiv.out$ROPE_high[[1]]
equiv.out %>%
  ggplot() +
  geom_hline(yintercept=rope_low, lty="dashed") +
  geom_hline(yintercept=rope_high, lty="dashed") +
  geom_pointrange(aes(x=Parameter, y=Hurricane, ymin=CI_low, ymax=CI_high, color=Estimate),
  position=position_dodge(.5)) +
  coord_flip() +
  scale_y_continuous(name="Hurricane Exposure", limits = c(-.3, .3), breaks=c(-.3,-.2,-.1,0,.1,.2,.3))+
  theme(legend.position="none")+
  theme_bw(base_size = 14)

ggsave(file="equivalency.png", path="Figures/", device='png', dpi=300)
